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ABSTRACT 



Aims. Fossil record methods based on spectral synthesis techniques have matured during the past decade, and their application to 
integrated galaxy spectra fostered substantial advances on the understanding of galaxies and their evolution. Yet, because of the lack 
^__^ of spatial resolution, these studies are limited to a global view, providing no information about the internal physics of galaxies. 

O Methods. Motivated by the CALIFA survey, which is gathering Integral Field Spectroscopy (IFS) over the full optical extent of 

600 galaxies, we have developed an end-to-end pipeline which: (i) partitions the observed data cube into Voronoi zones in order 
\^ to, when necessary and taking due account of correlated errors, increase the signal-to-noise ratio, (ii) extracts rest-framed spectra, 

including propagated errors and bad-pixel flags, (iii) feeds the spectra into the starlight spectral synthesis code, (iv) packs the results 
for all galaxy zones into a single FITS or HDF5 file, (v) performs a series of post-processing operations, including zone-to-pixel 
image reconstruction and unpacking the spectral and stellar population properties derived by starlight into multi-dimensional time, 
Q metallicity, and spatial coordinates. This paper provides an illustrated description of this whole pipeline and its many products. Using 

^ data for the nearby spiral NGC 2916 as a show case, we go through each of the steps involved, and present a series of ways of 

^ visualizing and analyzing this manifold. These include 2D maps of properties such as the velocity field, stellar extinction, mean 

^ ages and metallicities, mass surface densities, star formation rates on diff'erent time scales and normalized in diff'erent ways, plus ID 

I— ~i averages in the temporal and spatial dimensions, leading to evolutionary curves and radial profiles of physical properties. Projections 

of the stellar light and mass growth onto radius-age diagrams are introduced as a means of visualizing galaxy evolution in time and 
^"H space simultaneously, something which can also be achieved in 3D with snapshot cuts through the (x,y, t) cubes. 

^ Results. The results provide a vivid illustration of the richness of the combination of IFS data with spectral synthesis, of the insights 

OO on galaxy physics provided by the variety of diagnostics and semi-empirical constraints obtained, as well as a glimpse of what is to 

OO come from CALIFA and future IFS surveys. 

o 

|/-N Key words, galaxies: evolution - galaxies: stellar content - galaxies: general- techniques: imaging spectroscopy 

^^ 1 ■ Introduction galaxy. Surveys based on integrated spectra, oflTer a richer list of 

^-H diagnostics of the gaseous and stellar components, but lack spa- 

•^ Last decade technology has lead to massive surveys either in tial resolution and suffer from aperture effects (one spectrum per 

. 1^ imaging or spectroscopy modes that have been very successful in ^t^ject, centered at the nucleus and not covering the full galaxy). 

^ providing information of the spectral energy distribution or cen- Moreover, integrated galaxy spectra do not allow to, for exam- 

JH tral/integrated spectra for a large number of galaxies (e.g. 2dF, Pl^' isolate morphological components (bulge, disc, bars), map 

^ Folkes et al. 1999' SDSS York et al. 2000' COMBO- 17 Bell the effects of mergers, trace secular processes such as stellar mi- 

et al. 2004; ALHAMBRA, Moles et al. 2008; COSMOS, Ilbert, gration, and other features of galaxy formation and evolution 

et al 2009). These data have allowed significant insight on the which can only be observationally tackled with a combination 

integrated properties of galaxies, as well as on the evolution of of imagmg and spectroscopic capabilities, 
global quantities such as the mass assembly (Perez- Gonzalez et The near future will see a proliferation of Integral Field 

al. 2008) and the star formation history (SFH) of the universe Spectroscopy (IFS) surveys, which will allow a detailed look at 

(Lilly et al. 1996; Madau et al. 1996). However, these surveys physics within galaxies, as opposed to the global view offered 

are limited in either spectral or spatial resolution. Imaging pro- by integrated light surveys. CALIFA, the Calar Alto Legacy 

vides spatial information, but at the expense of coarse spectral Integral Field Area survey is a pioneer in this area (Sanchez 

coverage (typically broad band filters), limiting the amount of et al. 2012), and others will follow soon, like SAMI (Croom et 

information on the stellar populations, and providing little or no al 2012), VENGA (Blanc et al. 2009), and MaNGA (Bundy et 

information on the gas (emission lines) nor the kinematics of the al., in prep). CALIFA is obtaining IFS for 600 nearby galaxies 
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(0.005 < z < 0.03) of all types, mapping stellar populations, ion- 
ized gas and their kinematics across the full optical extent of the 
sources. The number of spectra to be generated is of the same 
order of the whole SDSS (~ 10^). 

Clearly, however, these will not be "just another million 
galaxy spectra". Processed through the machinery of fossil 
record methods, which recover the history of a stellar system 
out of the information encoded in its spectrum (Walcher et al. 
2011 and references therein), CALIFA will provide valuable 
and hitherto unavailable information on the spatially resolved 
SFH of galaxies. An example of what can be done was pre- 
sented by Perez et al. (2012), where we have used data on the 
first 105 galaxies observed by CALIFA in conjunction with the 
spectral synthesis code starlight (Cid Fernandes et al. 2005) to 
study the spatially resolved mass assembly history of galaxies. 
We find that, for massive galaxies, the inner parts grow faster 
than the outer ones, a clear signature of inside-out growth, and 
that the relative growth rate depends on the stellar mass of the 
galaxy, with a maximum relative growth eflftciency for interme- 
diate masses (~ 7 x lO^^M©). 

The application of spectral synthesis methods to IFS dat- 
acubes opens new ways of studying galaxies and their evolution. 
In practice, however, this apparently simple extension of work 
done for spatially integrated spectra involves a series technical 
details and methodological issues, as well as challenges in the 
visualization and representation of the results. 

The purpose of this paper is to present a detailed account of 
all steps in the processing of CALIFA (x, y, A) datacubes to the 
multidimensional products derived from the starlight analysis. 
We thus focus basically on methodology, leaving the exploration 
of results for subsequent studies. Paper II (Cid Fernandes et al. 
2013) presents a thorough study of the uncertainties involved 
in these products, including the eff'ects of random noise, spectral 
shape errors, and evolutionary synthesis models used in the anal- 
ysis. Despite the emphasis on CALIFA and starlight, the issues 
discussed here and in Paper II are of direct interest to spectral 
synthesis analysis of IFS data in general. 

The paper is organized as follows. Section |2] presents the 
data. Section [3] provides a detailed description of the pre- 
processing steps which prepare the reduced data for a full spec- 
tral fitting analysis. Section |4]reviews the starlight code and the 
ingredients used in our analysis. It also introduces the pycAsso 
pipeline, used to pack and analyze the results of the synthesis. 
Our main results are presented in Section [5] which uses the 
nearby spiral NGC 2916 (galaxy 277 in the CALIFA mother 
sample) as a show case. This long section starts with considera- 
tions on how to evaluate the quality of the spectral fits, and then 
proceeds to the presentation of the synthesis products in diff'er- 
ent ways, like 2D maps of the stellar light and mass, extinction, 
mean ages and metallicities, and star formation rates maps on 
diff'erent time-scales, ID maps on either the spatial or tempo- 
ral dimensions, and attempts to visualize the evolution of galax- 
ies in time and space simultaneously. The emphasis is always 
on methodological issues, but the examples naturally illustrate 
the richness of the multidimensional products of this analysis. 
Finally, our results are summarized in Section [6] 

2. Data 

The goals, observational strategy and overall properties of the 
survey were described in Sanchez et al. (2012). In summary, 
CALIFA's mother sample comprises 939 galaxies in the 0.005- 
0.03 redshift range, extracted from the SDSS imaging survey to 
span the full color-magnitude diagram down to M^ < -18 mag. 



and with diameters selected to match the field of the IFU in- 
strument (~ r in diameter). A thorough characterization of the 
sample will be presented in Walcher et al. (in prep). CALIFA 
will observe ~ 600 galaxies from this mother sample, applying 
purely visibility-based criteria at each observing night. The final 
sample will be representative of galaxies in the local universe, 
statistically significant and well selected. 

The data analyzed in this paper were reduced using the 
CALIFA Pipeline version 1.3c, described in the Data Release 
1 article (Husemann et al. 2013). Each galaxy is observed with 
PMAS (Roth et al. 2005, 2010) in the PPAK mode (Verheijen et 
al. 2004; Kelz et al. 2006) at the 3.5m telescope of Calar Alto 
Observatory and with two grating setups: the V500, covering 
from ~ 3700 to 7500 A at FWHM ~ 6 A resolution, and the 
V1200 (~ 3650-4600 A, FWHM ~ 2.3 A). The wide wave- 
length coverage of the V500 is ideal for our purposes, but the 
blue end of its spectra, which contains valuable stellar popula- 
tion tracers such as the 4000 A break and high order B aimer 
lines, is aff'ected by vignetting over a non-negligible part of the 
field of view. To circumvent this problem we worked with a com- 
bination of the two gratings, whereby the data for A < 4600 A 
comes from the VI 200 (which does not suff'er from this prob- 
lem at this wavelength range) and the rest from the V500. The 
spectra were homogenized to the resolution of the V500. 

For each galaxy the pipeline provides a (x, y. A) cube with a 
final sampling of (T', T', 2A). The cube comprises absolute cal- 
ibrated flux densities at each spaxel (F^), corrected by Galactic 
extincion using the Schlegel, Finkbeiner & Davies (1998) maps 
and the extinction law of Cardelli, Clayton & Mathis (1989) law 
with Ry = 3.1. Besides fluxes, the reduced cubes contain care- 
fully derived errors (6^), and a bad-pixel flag (/?i) which signals 
unreliable flux entries, due to, say, bad-columns, cosmic rays and 
other artefacts (see Husemann et al. 2013 for details). 

3. Pre-processing steps: From datacubes to 
STARLIGHT input 

A series of pre-processing steps are applied to the reduced dat- 
acubes, with the general goal of extracting good quality spectra 
for a stellar population analysis with starlight or any other sim- 
ilar code. These are: 

1 . Definition of a spatial mask 

2. Outer ^/A^ mask 

3. Refinements of the /?^ (flag) and 6^ (error) spectra 

4. Rest- framing 

5. Spatial binning 

6. Resampling in A 

This section describes each of these operations. 



3.1. Spatial masks 

The first step is to mask out foreground stars, artefacts and low 
S /N regions. As a starting point, candidate spurious sources are 
identified applying SExtractor (Berlin & Amouts 1996) to the 
SDSS r-band image. The masked image is then matched to the 
CALIFA field of view using the WCS information available in 
the headers to reproject the SDSS masks to the scale, orien- 
tation and pixel size of the CALIFA cubes. To improve upon 
this first approximation we perform a visual inspection of each 
galaxy, correcting its mask interactively whenever necessary. For 
instance, we manually fixed a few unmasked pixels which re- 
mained in the borders of each masked object (probably due to 
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Fig.l. (a) SDSS stamp for NGC 2916 (CALIFA 277). As in 
other images throughout this paper, the circles mark R = I and 2 
Half Light Radius (HLR). (b) CALIFA image in the rest-frame 
5635 ± 45 A interval, after applying the spatial mask, (c) Map 
of the S/N Sit 5635 A. (d) S/N map after the Voronoi binning. 

(e) Number of spaxels in each Voronoi zone. Note that no spatial 
binning (A^^ = 1) was necessary throughout the inner ~ 1 HLR. 

(f) Percentage of bad pixels in the 3800-6850 A range. North is 
up and east is towards the left in all images. 



a combination of the slightly better PSF of the SDSS images 
and the accuracy of the CALIFA WCS information used in the 
matching). Also, in a number of cases we found that the initial 
mask included bright Hii regions, while in others very dim point- 
like sources can be distinguished in the SDSS image, but have no 
clear counterpart in the CALIFA image nor spectra. In all cases, 
the masks were corrected by hand. 

The resulting masks are free from spurious spaxels. Unlike 
the other operations described below, this step could not be fully 
automated. For the size of CALIFA, this is still a feasible, albeit 
laborius, task. A last step in defining the spatial mask is to im- 
pose a spectral quality threshold based on the S/N of the data. 
This step has been fully automated in the pipeline described next. 

3.2. Flags, errors and segmentation: The qbick pipeline 

The next steps (items ii to vi) were all packed in a fully auto- 
mated package, named qbick, built specifically for CALIFA but 
sufficiently general to handle any FITS format to pre-process 
datacubes for a 3D spectral fitting analysis. 

The b;i flags produced by the reduction pipeline (/?^ = for 
good pixels and > for flagged ones) do a good job in iden- 



tifying problematic pixels. For uniformity, we enforce Z?^ > 
around the strongest sky lines seen in the CALIFA data: Hgl 
4358 A, Hgl 5461 A, 01 5577 A, Nal D (around 5890 A), 01- 
OH (6300 A), 01 6364 A, and the B-band atmospheric absorp- 
tion (Sanchez et al. 2007). We have further set /?i > when the 
value of the error is very large (specifically, when 6^ > 5F^), 
but this is merely a cosmetic choice, as these pixels would not 
influence the spectral fits anyway. Anomalously small 6^ values, 
on the other hand, can have a disproportionately large weight on 
the analysis. In the very few cases where this happens, the pixels 
were flagged. These refinements are useful to circumvent poten- 
tial problems in the spectral fitting analysis, specially when the 
whole data set is processed in "pipeline mode". 

After these fixes all spectra are rest-framed using the redshift 
derived by the reduction pipeline for the central 5'' spectrum 
(Sanchez et al. 2012), included in the original headers. This step 
includes a (1 -h z)^ correction for the energy, time and bandwidth 
changes with redshift. 

We used the (rest wavelength) window between 5590 and 
5680 A to estimate sl S/N ratio, defining the signal as the mean 
flux and the noise as the detrended standard deviatioqj of F^, 
both excluding flagged pixels [^Regions outside the outer S /N = 
5 iso-contour were added to the spatial mask, and discarded from 
the analysis. In a few cases this threshold was lowered due to low 
surface brightness. 

Fig.[l|) shows the 5635 ± 45 A image of CALIFA 277, after 
the final spatial mask. The hole above the nucleus corresponds 
to the foreground source seen in the SDSS image (panel a). The 
hexagonal pattern of the fiber bundle is clearly visible, since in 
this example S/N > 5 over most of the field of view (panel c). 

The 5635 A image is also used to compute the galaxy's Half 
Light Radius (HLR), adding up fluxes in spaxels sorted by dis- 
tance to the nucleus and picking the radius at which half of 
the total flux is reached. For CALIFA 277, HLR = 17.9''. This 
datacube-based HLR may not be the ideal measure of the size 
of a galaxy, but it is the most convenient metric for the present 
analysis for the simple reason that it is based on the exact same 
data. In Fig. [T] and all images throughout this paper circles are 
drawn at 7? = 1 and 2 HLR to guide the eye. 

The following step (item v in our list) consists of applying a 
segmentation structure to the data, grouping spaxels into spatial 
zones. All spaxels belonging to the same zone are assigned a 
common integer ID label, unique for each zone, qbick stores the 
zone-to-xj tensor needed to map zones to pixels and vice-versa. 
In this paper we use spatial binning to increase the S/N for the 
spectral analysis. Other applications may prefer to group spaxels 
according to, say, morphology or emission line oriented criteria. 
QBICK admits an externally provided segmentation map to allow 
for such choices. 

The zoning scheme adopted here is based on the Voronoi 
tesselation technique, as implemented by Cappellari & Copin 
(2003), but with modifications to account for correlated errors 



(Section 3.3 ). The code is feed with the signal and noise images 
in the 5635 ± 45 A range discussed above. The target S/N is 
set to 20. In practice, most individual spaxels inside 1 HLR sat- 
isfy this condition. Fig.[T]l shows the S/N map of CALIFA 277 
after the spatial binning. For this galaxy, the 3649 non-masked 
spaxels were grouped into 1638 zones, the overwhelming ma- 



^ The detrended std. dev. is the rms variation with respect to a linear 
fit, useful to remove the effect of the spectral slope in the window. 

^ Note that this is an apparently superfluous definition, since we al- 
ready have the noise spectrum. This "empirical ''S/N is nevertheless 



useful, as explained in Section 3.3 
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jority (1527) of which are not binned at all. As seen in the map 
of the number of spaxels per zone (A^^, Fig.[T^), spatial binning 
effectively starts at 7? ~ 1 HLR. Discounting the single spaxel 
zones, the median A^ is 10, but values up to 92 are reached in 
the outskirts. 

Once the spatial binning map is defined, the flux, error and 
flag spectra are computed for each zone. This operation must 
account for the flagged pixels. The total flux in zone z, containing 
A^ spaxels whose fluxes are F^,y^, is given by 



'A,z ^ 









(1) 



where 0^,y^ = for flagged pixels (b;i^k > 0) or 1 otherwise 
(i^A,k = 0), and O^^^ = 2 (pA,k is the number of useful spaxels 
at wavelength A in the zone. The zone flux is thus A^ times the 
average flux over all non-flagged pixels. When none of the F;i^k 
fluxes are flagged, this equation simply sums all the fluxes in 
the zone, whereas when this is not the case, flagged entries are 
replaced by the average of non-flagged ones. 

This scaling recipe is reasonable as long as the number of 
flagged pixels is not comparable to A^. For example, in a zone 
of A^ = 10 spaxels, one would not like to take the resulting F^,^ 
seriously if 8 of them are flagged. To deal with this issue we flag 
away entries whenever less than half of the spaxels in a zone 
have credible data at the given A: 



n,z 



ifO^,,>A,/2; 

1 ifO^,, <A,/2. 



(2) 



For our example galaxy, the fraction of flagged pixels in the 
3800-6850 A interval used in the spectral synthesis analysis 
ranges from 5 to 8%, with a median of 6% (see Fig.[T]0. 

The computation of the error in F^^^ follows a recipe similar 
to eq. [T] but adding in quadrature 



N ^' 



(3) 



where y^^ accounts for correlated errors (see below). 

QBicK also produces the total spatially integrated spectrum, 
following these same prescriptions. This is useful to compare 
the analysis of the whol e to the sum of the analysis of the parts 
(which we do in Section |5]9]). 

Finally, all the spectra are resampled to 2 A. The resulting 
zone spectra are stored in ASCII and/or FITS files, with flux , 
error and flag columns necessary for a careful spectral analy- 
sis. All the quality control images produced in this process are 
packed in a single FITS file along with the technical parameters 
adopted in the different steps (S/N clipping, segmentation map, 
etc.) 

3.3. Spatial binning and correlated errors 

As already mentioned, the data cubes have a spatial sampling 
of r' X r'. Because of the three-fold dithering pattern used in 
the observations, each spaxel has a varying contribution from a 
number of fibers. Since one fiber can contribute to more than one 
spaxel, the noise in a given spaxel is partially correlated with that 
of its neighbors. This implies that when spectra from adjacent 
spaxels are added together the noise in the resulting spectrum 
has to be calculated taking spatial covariances into account. As 
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Fig. 2. Ratio of the measured over the idealized noise for Voronoi 
zones comprising A^ spaxels. The idealized noise is that derived 
under the assumption of uncorrelated errors, while the measured 
one is obtained from the detrended standard deviation of fluxes 
in a ±45 A range around 5635 A. Gray dots represent yS^ values 
for 99267 spaxels from 109 galaxies. Circles show the median yS^ 
for each A^, and the solid line shows equation]?] The top panel 
show the histogram of A^-values (notice the logarithmic scale). 



an illustrative extreme example, consider A^ completely corre- 
lated (ie., ~ identical) spaxels in a zone. Clearly, co-adding their 
spectra would not result in the desired increase in the S/N, as 
both signal and noise scale linearly with A^. Unless told other- 
wise, however, the zoning algorithm assumes that the spaxels 
are independent, and hence that the noise in the zone is given by 



V 



2 6^, resulting in an improvement of the S/N of the co-added 

spectrum by ~ VA^, whereas in fact it has not changed at all|j 

An empirical correction was devised to account for the 
global behaviour of the spatial correlation of the signal in 
CALIFA data. The way we proceed is as follows: First we define 
initial Voronoi zones with the Cappellari & Copin (2003) algo- 
rithm. Because it assumes uncorrelated input, the code presumes 

the zone error is J^j ^l- We then measure the "real" noise {e^) in 
the resulting zone summed spectra from the (detrended) rms in 
the 5635 ±45 A range. The ratio of these two numbers, which we 
cdiWPz, is represented versus the number of spaxels in each zone 
in Fig. [2] Small grey dots arcyS^ values for the individual zones, 
circles represent the median value computed at each A^, while 
the shaded region represents the 16 to 84 percentiles. There are 
nearly 10^ points in this plot, obtained from 244196 spaxels in 
109 galaxies (all those observed with the V500 and VI 200 se- 
tups up to May 2012). The top panel shows the histogram of A^. 



^ 6yt denotes the noise in spaxel k, defined as the detrended rms in the 
same 5635 + 45 A range used to define the signal for spatial binning 
purposes. 
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The rise and flattening of the Pz-^z relation reflects the ex- 
pected behavior. Zones with small A^^ group adjacent and highly 
correlated spaxels, producing the initial rise of p^. As farther 
apart spaxels are aggregated, the errors become less correlated, 
leading to the flattening observed. The scatter in the relation is 
likely related to the way the Voronoi algorithm assemble spaxels. 
Also, the three-fold dithering pattern is not completely uniform, 
so the amount of correlation is not exactly the same across the 
face of the fiber bundle. Still, theyS^(A^^) relation portrayed in Fig. 
|2]is tight enough to derive an excellent statistical correction. 

The solid line represents the fit of the following function: 



PziNz) = 



15A^, 



15-hA^, -1 



1/2 



(4) 



where the fit was carried out considering only the filled circles. 
We can now apply this correction directly in the zoning code, 
replacing its idealized computation of the noise by 



ez=Pz(Nz) 



\ 



Z 



ek' 



(5) 



which accounts for the eff'ects of spatially correlated errors. This 
correction was applied /I by /I (cf. equation [Sj. An indepen- 
dently derived correction very similar to equation]?] is presented 
in Husemann et al. (2013). 

This scheme naturally leads to larger zones than those ob- 
tained under the (erroneous) hypothesis of independent spaxels, 
but ensures that the target S jN of 20 is actually achieved in all 
but the outermost zones (when the algorithm runs out of spaxels 
to bin). Typically, we obtain 1000 zones per galaxy. 



4. STARLIGHT runs and PyCASSO 

STARLIGHT is a spcctral synthesis method which combines the 
spectra from a base in order to match an observed spectrum O^- 
The search for an optimal model M^ and the coefficients of the 
population vector x associated with the A^^ base elements also 
allows for reddening, a velocity off-set v^ and a stellar velocity 
dispersion cr^. The code was first introduced in Cid Fernandes 
et al. (2004) and substantially optimized by Cid Fernandes et al. 
(2005). The code, didactic introductions and a user manual are 
available at www.starlight.ufsc.br. We are actually using a new 
and yet to be documented version of starlight, but since most of 
the new features are not really used in our analysis, the publicly 
available manual remains an accurate reference for the purposes 
of this paper. 

The virtues and caveats of full spectral fits are discussed in 
a number of articles (e.g., Panter et al. 2003; Ocvirk et al. 2006; 
Cid Fernandes 2007; Tojeiro et al. 2007; Mc Arthur, Gonzalez & 
Courteau 2009; Sanchez-Blazquez et al. 2011). A general con- 
sensus in the field is that uncertainties in the results for indi- 
vidual objects average out for large statistical samples (Panter 
et al. 2007). With CALIFA, each galaxy is a statistical sample 
per se ! Hence, even if the results for single spaxels (or zones) 
are not iron-clad, the overal trends should be robust (see Paper 
II for a detailed evaluation of the uncertainties in our starlight- 
based analysis of CALIFA data.) It is also worth pointing out 
that despite the diversity of spectral synthesis methods, substan- 
tial changes on the results are more likely to come from revisions 
in the input data (as happened, for instance, with the recalibra- 
tion of SDSS spectra between data releases 5 and 6) and from 



updates in the base models, the single most important ingredient 
in any spectral synthesis analysis. 

The results reported in this paper rely on a spectral base 
of Simple Stellar Populations (SSPs) comprising 4 metallici- 
ties, Z = 0.2, 0.4, 1 and 1,57©, and 39 ages between t = 10^ 
and 1.4 X 10^^ yr. This base (dubbed ''GAf' in Paper II) com- 
bines the Granada models of Gonzalez Delgado et al. (2005) for 
t < 63 Myr with those of Vazdekis et al. (2010, as updated by 
Falcon-Barroso et al. 2011) for larger ages. They are based on 
the Salpeter Initial Mass Function and the evolutionary tracks 
by Girardi et al. (2000), except for the youngest ages (< 3 Myr), 
which are based on Geneva tracks (Schaller et al. 1992; Schaerer 
et al. 1993a,b; Charbonnel et al. 1993). A comparison of results 
obtained with other bases is presented in Paper II, to which the 
reader is also referred for examples of the spectral fits. 

The fits were carried out in the 3800-6850 A interval. 
Reddening was modeled with the Cardelli, Clayton & Mathis 
(1989) curve, wit 7?y = 3.1. As in Cid Fernandes et al. (2005), 
the main emission lines as well as the Nal D doublet (because 
of its sensitivity to ISM absorption) were masked. All runs were 
performed in the lAA-GRID, a network of computers which al- 
lows us to process 10^ spectra in less than a day. 

As documented in the user manual, starlight outputs a large 
array of quantities. In broad terms, these can be categorized as 

1. Input data: File names, configuration options, and other in- 
formation provided either by the user explicitly or derived 
from the user-provided spectrum. Base-related data are also 
reported, including ages and metallicities of each compo- 
nent, as well as the corresponding light-to-mass ratios (at the 
chosen normalization wavelength. An) and a retumed-mass 
correction factor. _ 

2. Fit results: Figures of merit (x^, A), kinematic parameters 
(v*, (T*), the V-band extinction (Ay), total stellar masses and 
luminosities, and population vectors, expressed in terms of 
light (x) and mass (fi) fractions. 

3. Spectral data: The input (O^), best model fit (M^), and 
weight spectra (w^, which equals e^^ except for ffagged, 
masked and clipped pixels). 

All these data are stored in a plain ASCII file, one for each 
fitted spectrum. The population vectors can be post-processed in 
a number of ways to obtain SFHs, to compute the mass growth 
as a function of time, to condense the results into first moments 
like the mean t and Z, etc. This scheme works fine for studies of 
independent objects, like SDSS galaxies or globular clusters. 

IFS data, however, require a more structured level of orga- 
nization, as it is cumbersome and inefficient to handle so many 
files separately. To handle datacubes we developed pycAsso, the 
Python Califa Starlight Synthesis Organizer. pycAsso comprises 
three main parts: 

1 . A writer module, which packs the output of all starlight fits 
of the individual zones of a same galaxy into a single FITS 
(or HDF5) file, which also contains all information gener- 
ated during the pre-processing steps (Section [3]), as well as 
information propagated from the original datacube. 

2. A reader module, which reads this file and structures all the 
data in an easy to access and manipulate format. 

3. A post-processing module, which performs a series of com- 
mon operations such as mapping any property from zones 
to pixels, resampling and smoothing the population vectors, 
computing growth functions, averaging in spatial, time or 
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Fig. 3. Maps of spectral-fit quality indicators, (a) Mean relative 
deviation A, (b) x'^ P^r fitted flux, (c) number of fluxes used in 
the fit, and (d) number of clipped fluxes. 

metallicity dimensions, etc. All these functionalities are con- 
veniently wrapped as python packages, easily imported into 
the user's own analysis code. 

pycAsso will be the main work-horse behind a series of arti- 
cles by our collaboration, so we dedicate the remainder of this 
paper to illustrating its main products. 

5. Results 

The Sb galaxy NGC 2916 (CALIFA 277) was chosen as a show- 
case. For the adopted distance of 56 Mpc, V (= 1 CALIFA 
spaxel) corresponds to 270 pc. With an r-band absolute mag- 
nitude of -2 LI 7, w - r = 2.18 and g - r = 0.54 colors, and 
concentration index C = 2.16, NGC 2916 can be characterized 
as a blue-disk galaxy with an intermediate mass (~ 10^ ^Mo). 
Rudnick, Rix & Kennicutt (2000) have shown that the galaxy 
show some degree of lopsidedness, indicating a weak interaction 
(possibly with its irregular satellite 5' away, Gutierrez, Azzaro 
& Prada 2002) or a minor merger. Moorthy & Holtzman (2006) 
have obtained long-slit spectra and measured emission lines and 
Lick indices along the position angle PA = -80°. They found 
that the nuclear spectrum is well within the AGN wing of the 
BPT diagram [Oiii]/Uj3 vs. [NiiJ/Hof, in a location close to what 
is now accepted as borderline between Seyferts and LINERs 
(Kewley et al. 2006). 

Our Voronoi binning of the CALIFA datacube produces 
N = 1638 zones. The corresponding N -\- I spectra (the -hi cor- 
responding to the spatially integrated spectrum) were fitted with 
STARLIGHT and the results feed into pycAsso. The following sub- 
sections describe how we handled the data and the products of 
the analysis. 

5.1. Fit quality assessment 

As with ID spectra, the quality of the spectral fits should be ex- 
amined before proceeding to interpretation of the results. This 



subsection describes some standard methods to do this, the 
caveats involved and strategies to circumvent them. 
Fig.|3^ shows the map of the fit quality indicator 



vl^ 



1 



Nfr 



M, 



M, 



(6) 



where O^ and M^ are the observed and model spectra, respec- 
tively, and the sum is carried over the A^^^ wavelengths actu- 
ally used in the fit, ie, discarding masked, flagged and clipped 
pixels |j For this galaxy A spans the range between 1.0 and 
5.6%, with a median value of 3.0%. A increases towards the 
outer regions, where the S IN is smaller (Fig. [ip). The x^ = 
2 w^(O^-M^)^ map (Fig.pp) shows the opposite behavior, since 
the errors also increase outwards. Because of its non-explicit re- 
liance on the (often hard to compute, if not altogether unavail- 
able) 6^ spectrum, as well as because of its easily grasped mean- 
ing, A is a more useful figure of merit to assess fit quality. 

Inspection of the highest A spectra often reveal non-masked 
emission lines or artifacts, like imperfect masking of foreground 
sources, slightly misplaced b^ flags, or (9^ values which should 
have been clipped by starlight but were not because of large 6^. 
We emphasize that such bugs are very rare in CALIFA. The me- 
dian A for the -10^ zone spectra analysed so far is just 4% (cor- 
responding to an equivalent S/N of 25), and in less than 2% of 
the cases A exceeds 10%. This high rate of success is mainly due 
to the carefully derived errors and flags in the reduction pipeline 
and further refined in our pre-processing steps (Section [3]). 

One can use such maps to define spatial masks over which 
the spectral fits satisfy some quality threshold. Experiments with 
the whole data set suggest that a reasonable quality-control limit 
should be in the A = 8 to 10% range. Not surprisingly, the large 
residuals tend to be located in the outer regions of the galaxies. 
Generally speaking, results beyond 2-3 HLR should be inter- 
preted with greater care. No A-based cut is needed for our exam- 
ple galaxy. For the benefit of starlight users, it is nevertheless 
worth making some further remarks on quality control. 

5.1 .1 . Subtleties and caveats on quality control 

In full spectral fitting methods, seemingly trivial operations like 
imposing a fit-quality threshold often hide subtleties which can 
easily go unnoticed, particularly when processing tons of data. 
The paragraphs below (focused on starlight but extensible to 
other codes) review some of these. 

First, one should distinguish cases where A is large because 
of bad data from those where large residuals arise because of the 
user's failure to properly inform, through spectral masks (m^) or 
flags (b;i), spectral regions which should be ignored in the fit. 
For instance, one can have an excellent spectrum of an HII re- 
gion with several strong emission lines besides those included in 
a generic emission line mask feed into starlight, causing an ar- 
tificially large A. Large velocity off'sets can have a similar eff'ect, 
as emission lines get shifted out of fixed m^ windows. Similarly, 
sky residuals and other artefacts missed out by the /?^ flags lead 



"^ Notice that eq.[6]is nearly, but not exactly identical to what is called 
"adev" in starlight's output, the difference being the use of M^ instead 
of Oa in the denominator. This avoids faulty (but not flagged) O^ en- 
tries with very low flux aff'ecting the statistics. This rarely occurs in 
CALIFA, and when it does the pixels invariably have large errors, and 
so are irrelevant for the fit, but can still aff'ect the value of A. 
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Fig. 4. (a) Synthesitc surface brightness at 5635 A, (b) stellar 
extinction map (Ay), (c) dereddened Xi5635 image, (d) stellar 
mass surface density. 



to large A even when the overall spectrum is good. In short, not 
everything which fails a blind quality control is actually bad. 

The clipping options implemented in starlight help spotting 
pixels which are too hard to fit and thus probably represent non- 
stellar or spurious features. Our fits use the "NSIGMA" clipping 
method and a conservative Act threshold, meaning that we only 
clip pixels when |(9^ - M^| > 46^. Fig. [3]: maps the number of 
clipped points in CALIFA 277. One sees a peak in the central 
pixels, where the S IN is so high that our Act clipping does not 
forgive even small O^ - M^ deviations (most often associated 
with problems with the base models rather than with the data). 
Elsewhere, very few pixels were clipped. The point to highlight 
here is that clipping only worked because we have a reliable e^. 
Other clipping methods can be used when this is not the case, but 
in all cases the user should always check carefully the output, as 
it may easily happen that too many points are clipped. 

Custom-made spectral masks also help. For instance, since 
Mateus et al. (2006), the starlight analysis of SDSS spectra 
employs individual masks constructed by searching for emission 
lines in the Ox - Mx residual spectrum obtained from a first fit, 
and taking into consideration the local noise level. Ox is then 
refitted with this tailor-made m^, circumventing some of the is- 
sues above. This refinement has not yet been implemented for 
CALIFA. 

Errors, flags and masks are, of course, secondary actors in 
any spectral synthesis analysis, but these general remarks illus- 
trate that it pays ofl' to dedicate them some attention. 



5.2. 2D maps: Stellar light, mass and extinction 

One of the main products of spectral synthesis is to convert light 
to mass. Fig.|4]illustrates the results for CALIFA 277. Its top-left 
panel shows the surface density of the luminosity at A^ = 5635 
A, while the top-right panel show the derived extinction map. 
The dust corrected image and the mass surface density are also 
shown. 



The efl'ects of the Voronoi binning are present in all panels, 
but are much more salient in the Ay map. This is because the 
light and mass images were ''dezonified" by scaling the value at 
each xy spaxel by its fractional contribution to the total flux in 
a zone (z). For instance, for the mass surface density (At) this 
operation reads 



Atx 



Axy A^ 



XWr 



(7) 



where Axy (A^) denotes the area in a spaxel (zone), and 



Zjxylz^ xy 



(8) 



with Fxy as the mean flux in the 5635 ±45 A region, the same 
used in the Voronoi zoning (Section [3]). This operation was ap- 
plied to luminosity and mass related quantities, producing some- 
what smoother images than obtained with Wxyz = 1. Intensive 
properties (like Ay, mean ages, (Ti, and etc.), however, cannot be 
dezonified. 

The stellar extinction map shows low values of Ay (of order 
0.1-0.2 mag), with slight enhancements in the nuclear region 
and the arms. Overall, however, there is relatively little variation 
across the face of the galaxy (see also Fig. [TO]). The light and 
mass images have a similar structure, both showing the bulge at 
R ^ 0.5 HLR, and the disc beyond that. The spiral arms are less 
prominent in terms of mass than in light because of the higher 
L/M of stars in the arms. 

The total stellar mass obtained from the sum of the zone 
masses is M = 6.5 x lO^^M©. This is the mass locked in stars 
nowadays. Counting also the mass returned by stars to the in- 
terstellar medium, M' = 9.0 x lO^^M© were involved in star 
formation.These values ignore the mass in the masked region 
around the foreground star northeast of the nucleus. pycAsso can 
fill in such holes with values estimated from (circular or ellip- 
tical) radial profiles. For CALIFA 277, this correction increases 
the masses quote above by 4%. 



5.3. 2D maps: Kinematics 

Fig. [5] shows the v* and cr^ fields estimated by starlight, as 
well as a position- velocity diagram and a radially averaged cr^ 
profile. The v^ field indicates a projected rotation velocity of 
~ 200 km/s. The V500 resolution prevents the determination of 
(T* values below -150 km/s, and the flat cr^(R) profile beyond 
~ 0.7 HLR reflects this resolution limit (notice that we have not 
corrected the values for the instrumental resolution in this plot). 
Only in the inner regions the cr^ values are trustable, reaching 
200 km/s in the nucleus (equivalent to 140 km/s after correcting 
for instrumental broadening). 

The kinematical information derived from our analysis will 
be superseded by studies based on the higher spectral resolution 
VI 200 datacubes (Falcon-Barroso et al, in prep). Eventually, one 
can envisage feeding the parameters derived from these more 
precise analysis back into the starlight fits (using its fixed- 
kinematics mode). In fact, this feedback may well turn out to 
improve our stellar population analysis, given the potential de- 
generacy between cr^ and Z (decreasing the former while in- 
creasing the latter have the same global eff'ect of making metal 
lines deeper; Koleva et al 2009; Sanchez-Blazquez et al. 2011). 
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Fig. 5. Kinematical products of the spectral fitting. The cr^ val- 
ues are not corrected for instrumental broadening, which dom- 
inates below 140 km/s. The horizontal stripes in the position- 
velocity and cr^(R) panels correspond to Voronoi zones. 



5.4. 2D maps: Mean ages and metallicities 

The crudest way to quantify the SFH of a system is to compress 
the age and metallicity distributions encoded in the population 
vectors to their first moments. For this purpose we will use the 
following definitions: 



(log Ol = ^ Xtz log t 
t,z 

t,Z 



(9) 



(10) 



where Xtz is the fraction of light at A^ due to the base population 
with age t and metallicity Z. The mass weighted versions of these 
indicators, (log Om and (Z)m, are obtained replacing Xtz by the 
mass fraction yu^z- 

Fig. [6] shows the light and mass weighted mean (log) age 
and metallicity maps. The (log t}i image shows a steady increase 
towards the center. Outside 1 HLR, traces of the spiral arms are 
noticeable as regions of lower age (as in the SDSS color image 
in Fig. [T^, the arms are brighter the western half of the image). 
Because of the large weight of old populations, (log Om spans 
a smaller dynamical range than (log Ol, hence producing lower 
contrast maps, but the outwardly decreasing age is still visible. 
Negative gradients are also clearly present in metallicity, with 
indications of flattening within the bulge region. As in other 2D 




Fig. 6. Luminosity (top) and mass (bottom) weighted mean age 
(left) and metallicity (right) maps for CALIFA 277. 



maps, the small scale fluctuations towards the edges are at least 
in part due to the lower S/N. 

The fact that t and Z grow in the same sense suggests that 
the results are not badly afl'ected by the age-metallicity degener- 
acy. As shown by Sanchez-Blazquez et al. (2011), full spectral 
synthesis is much less sensitive to this than conventional index- 
based approaches. Fig. [7] plots extinction versus mean age ver- 
sus metallicity, the main properties involved in spectral synthe- 
sis. Points are color and symbol coded by their distance from 
the nucleus. One sees that the highest values of Ay are found 
for old central (R < 0.5 HLR) and young outer (R > 1 HLR) 
populations, and that Ay bears no correlation with Z. The mid- 
dle panel shows the positive t-Z relation inferred from the 2D 
maps. Within 0.5 HLR (green diamonds), however, t anticorre- 
lates with Z, most likely due to the age-metallicity degeneracy. 
The mean Z values in this central region straddle 1 and 1.5 Z©, 
the two highest metallicities in our base. 

5.5. 2D maps: xy, xi, xq 

Population synthesis studies in the past found that a useful way 
to summarize the SFH is to condense the age distribution en- 
coded in the population vector into age ranges. This strategy 
comes from a time when the analysis was based on equivalent 
widths and colors (Bica 1988; Bica et al. 1994; Cid Fernandes 
et al 2001, 2003), but it was also applied to full spectral fits 
(Gonzalez Delgado et al. 2004). 

Fig. [8] presents images of the light fraction due to Young, 
Intermediate and Old populations (xy, xj and xq), defined as 
those with t < 0.14, 0.14 < ^ < 1.4, and t > 1.4 Gyr, respec- 
tively. As usual, the choice of borderlines is somewhat subjec- 
tive, constrained only by the underlying idea of grouping base 
elements covering relatively wide age ranges. For the base used 
in our starlight runs the Y, I and O bins contain 4 x 20, 4x10 
and 4x9 SSPs, respectively (where the 4x comes from the 4 
metallicities). 

The plots show that the spiral arms stand out more clearly in 
the Xy map, specially in the western half of the galaxy, as also 
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Fig. 7. Mean age versus mean metallicity versus extinction, cod- 
ing zones by their distance to the nucleus. Green diamonds: 
R < 0.5 HLR. Red triangles: 0.5-1 HLR. Blue stars: 1-2 HLR. 
Black dots: > 2 HLR. 



seen in the SDSS color composite (Fig. [T^). Very little of the 
light from R < 1 HLR comes from young stars. Intermediate 
age populations do contribute more, but the central ~ 0.5 HLR 
is completely dominated by old stars. Overall, however, despite 
the added informational content, these maps do not visibly add 
much to the (log Ol image (Fig.|6]). 
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5.6. 2D maps: Star formation rates and IFS-based variations at lookback time ^PlEq 
over Scaio's b parameter 



Fig. 8. Maps of the percentage contribution of Young (t < 0.14 
Gyr), Intermediate age (0.14-1.4 Gyr), and Old (> 1.4 Gyr) stars 
to the observed light at 5635 A. 



where Mf-^y is the mass (at xy and per unit area) of stars formed 



11 



The base used in our fits comprises instantaneous bursts (ie., 
SSPs), so, despite the large number of ages considered, our SFHs 
are not continuous and hence not derivable. There are, nonethe- 
less, ways of defining star formation rates (SFR). 

The simplest way to estimate a SFR is to cumulate all the 
stellar mass formed since a lookback time of tsF and perform a 
mass-over-time average: 



gives the mean SFR surface density 
since t = tsr- One can tune tsr to reach diff'erent depths in the 
past, but as tsr increases this estimator becomes increasingly 
useless, converging to the mass density divided by ^oo = 14 Gyr 
(the largest age in the base). 

It is often more useful to consider SFRs in relation to some 
fiducial value, instead of absolute units. The classical example 
is Scaio's birthrate parameter, b, which measures the SFR in 



tsF f-r 



(11) 



KtsF 



^ The stellar mass in this equation differs from that in eq. [t] Fig. 
[4] and the definitions of (log^M and {Z)m- The latter are corrected by 
the mass returned to the medium by stars during their evolution (thus 
reflecting the mass currently in stars), while in the computation of SFR's 
this correction should not be applied. We distinguish the two types of 
stellar masses with the prime superscript (Al versus At')- 
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Fig. 9. Spatially resolved SFR surface density (eq. 1 1 ) on the last 
tsF = 14 (top), 142 (middle) and 1420 (bottom) Myr, in units of 
different reference values. The left panels compare STUxyitsr) 
to ST^xyitoo), i.e, the all-times average at xy, thus providing 
a local version of Scalo's b parameter (eq. [12]). Middle panels 
compares the local STU to the all-times average over the whole 



galaxy (eq. 13 ). Panels on the right compare STIi to that of the 
galaxy as a whole over the same time scale (eq.[T4]). All images 
are on log scale from log 1 to log 5, such that only STUxy^tsr) 
values above the corresponding reference value are visible. 



the recent past {t < tsr) with respect to its average over the 
whole lifetime of the system|jThis is trivially obtained dividing 
STUxyitsp) by its asymptotic value 



h^'xyitsF) = 



SrKydsF) too I.t<ts,M',,y 



SrHxyitoo) tsF Mxy 



(12) 



The superscript "loc" is to emphasize that the reference lifetime 
average SFR is that of same spatial location xy. One should re- 
call that while young stars in a spaxel were probably born there 
(or near it), old ones may have migrated from different regions. 
Thus, despite the identical definitions, the physical meaning of 
b^^y is not really the same as for galaxies as a whole. 

IFS data allow for other definitions of b. For instance, one 
might prefer to compare STI^xyitsr) to the ST1^{t oo) of the 
galaxy as a whole, the bulge, the disc or any general region. This 
variation measures the "present" (t < tsr) "here" (at xy) to the 
"past" (^ < ^oo) in a spatial region 1i\ 



b'i.itsF) 



STKyitSf) 



STKniU 



(13) 



A formally similar but conceptually different definition is ob- 
tained using for the reference value in the denominator the SFR 



^ A related index is the so called specific SFR, which differs from b 
just by the ^oo factor in equation H^ 



surface density within the same time-scale but evaluated in a dif- 
ferent region: 



C%{tSF) = 



srn.yitsF) 



STKnitsf) 



(14) 



i.e., to compare not "now versus then", but "now here versus 



now there". Together with equations 12 and 13 this definition 
cover the different combinations of time and space enabled by 
the application of fossil methods to IFS data. 

Fig. [9] shows maps of b^^^ (left panels), Z?^ (middle), and c^ 
(right) for three values of tsr'- 14 (top), 142 (middle) and 1420 
Myr (bottom). In all panels the color scale is deliberately satu- 
rated to highlight regions where ST^itsp) is larger then chosen 
reference value — the dynamical range of the images goes from 
1 to 5 (0 to 0.7 in log) in the corresponding relative units. 

Panel a shows that spaxels which have formed stars over the 
past 14 Myr at a larger rate than their respective past average 
are all located in the outer regions. Nowhere within 7? ^ 0.5 
HLR one finds Z?^''''(14Myr) > 1. Considering the past 142 Myr 
(panel d), the inner "deficit" covers as much as 1.2 HLR. On the 
longest time scale considered (g) one again sees that the inner 
spaxels have been less active than their life-long average. 

Maps look considerably different in the middle column of 
panels, where the local STU is normalized to the past aver- 
age over the whole galaxy (eq.[T3]). The spiral arms of CALIFA 
277 appear better delineated in Z7^^^(14Myr) (panel b) than in 
any b^^^ map. Interestingly, the b^^^ map is practically featureless 
for tsF = 142 Myr (e). This suggests that star- formation is not 
continuous over this time scale, but instead happens in a bursty 
mode. Over the past 1.4 Gyr (h), one again sees traces of the 
arms, but with a lower amplitude than in panel b, qualitatively 
consistent with the cumulative effect of an intermittent sequence 
of short duration bursts. (Due to the logarithmic age resolution 
of fossil methods, such short bursts can only be recognized as 
such at the very young ages sampled in panel b.) 

The right column panels show our novel relative SFR in- 
dex c (eq. W^, using the whole galaxy as reference region. 



Unlike Scalo's b, this index does not compare present to past, 
but present to present elsewhere. Its 2D maps can be understood 
as "snapshots" of the SFR in the galaxy taken with an "expo- 
sure time" tsF' Keeping the "diaphragm" open for only the last 
14 Myr (panel c) highlights the ongoing star- formation in the 
galaxy. The spiral arms are clearly visible, and the structures are 
more focused than in the M^^(14Myr) map (panel b). Integrating 
for 142 Myr (f), the inner parts of the arms fade, but their outer 
(R > I HLR) parts brighten up in terms of c^^K Over the past 
142 Myr, these regions have formed more stars than anywhere 
else in the galaxy, even though comparing with panel e we find 
that only parts of the c^^^ > 1 regions are forming stars at a rate 
per unit area larger than that of the galaxy as a whole through 
its entire life (ie, b^^^ > 1). For tsF = 1-4 Gyr (panel i) the im- 
age becomes more blurred. Because of the long time scale, the 
c^^^(1.4Gyr) map is nearly indistinguishable from that obtained 
with b^^^ for the same tsF (panel h). 

One thus sees that, despite some degree of redundancy (the 
numerator is always the same), these definitions offer different 
and complementary views of the star formation in a galaxy. 

5.7. 1D spatial maps: radial profiles 

The information contained in 2D maps like the ones shown 
above can sometimes be hard to absorb. Azimuthal averaging 
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Fig. 10. Radial profiles of several properties. Grey points correspond to values in individual xy spaxels. The green circles and 
solid lines represent the average of the plotted quantity, while the error bars represent the dispersion in the 7?-bin. The dashed 
magenta horizontal line (only in the plots not involving surface densities) marks the value derived from the starlight analysis of the 
integrated spectrum, i.e. collapsing the xy dimensions of the datacube. The stars in the bottom panels mark the 4 metallicities in the 
base models. 



is a useful way to compress and help digesting 2D data. pycAsso 
provides both circular and elliptical xy-to-R conversion tensors. 
The examples below are based on a circular mapping. 

For quantities like At, one can think of two types of radial 
averaging. The first is to add up all the stellar mass in xy spaxels 
within a given 7?-bin and then divide by the bin area (counting 
only non-empty spaxels). The second is to average the surface 
density values for each spaxel directly. The same applies to, for 
instance, (log Ol (eq.[9]): One can either add up the value of prod- 



uct Ltxy X log t for all spaxels at a given R and then divide by the 
corresponding Y^Lxy, or else simply average the (log^L values 
for each spaxel in the 7?-bin. The first method, which we may call 
area averaging, eff'ectively collapses the galaxy to ID, but this 
kind of averaging cannot be applied to quantities which do not 
involve light or mass, like Ay and cr^. For the sake of uniformity, 
in this subsection we chose the second type of radial averaging 
(henceforth spaxel averaging), noting that the two methods al- 
most always give nearly identical results. 
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Fig. [T0| shows several pycAsso products as a function of R. 
The grey points represent values for individual spaxels. The solid 
line and green circles show the mean value of the plotted quan- 
tity among the spaxels in the same R bin, and the error bars map 
the corresponding dispersion. (As discussed in Paper II, actual 
error bars in radial profiles are much smaller due to the large 
number statistics.) Several of the features discussed while de- 
scribing the 2D images are clearly seen in these plots. The (log t) 
and (Z) gradients, in particular, are cleanly depicted in these 
plots. As pointed out in Section 5A_ means ages and metallic- 
ities increase ~ simultaneously towards the nucleus, while Ay 
does not change much, indicating that the classical degenera- 
cies in stellar population studies are not playing a strong role in 
shaping the results, with the possible exception of the central re- 
gions. As in Fig. [5 horizontal stripes of gray points in panels c, 
e, f, g and h of Fig. 10 correspond to spaxels in the same Voronoi 
zones. The stripes disappear in panels a, b and d because these 
extensive quantities were "dezonified" (cf. eq.jT]). 

5.8. 1D temporal maps: evolutionary curves 

The main power of fossil record methods is to open the time do- 
main, allowing evolutionary studies. The caveat is that it does so 
with the typical logarithmic resolution inherent to stellar evolu- 
tion, but several studies show that much can be learned about 
galaxy evolution despite this fundamental constraint. Clearly, 
though, our 39 ages base is highly overdimensioned, so some 
sort of SFH compression is necessary. As discussed by Asari et 
al. (2007), age resolutions between 0.5 and 1 dex are reasonable. 
We thus smooth the light and mass population vectors with a 
gaussian of FWHM = 0.5 dex in log t. 

Fig.[TT]plots a number of our synthesis products against age. 
To keep an at least partial representation of the spatial informa- 
tion, we plot evolutionary curves derived for diff'erent radial re- 
gions: The nucleus, defined as the central pixel (plotted in grey 
dashed lines), R < 0.5 (dashed green), < 1 (solid red), and > 1 
HLR (solid blue). 

The top panels show the growth of luminosity (Fig.[TT^) and 
mass (b)Q Among other things, this plot illustrates the diff'er- 
ence between light and mass, in the sense that regions inside and 
outside R = I HLR do not have the same mass, despite hav- 
ing (by definition) the same light. While HLR =4.8 kpc, the 
radius containing half of the mass is 3.4 kpc. This is a direct 
consequence of the stellar population gradients in this galaxy. 
Gonzalez Delgado et al. (in prep) investigates the relation be- 
tween light and mass based radii for diff'erent types of galaxies. 

It is also seen that both light and mass grow at different 
speeds for different regions. This is better appreciated in pan- 
els c and d, were each growth curve is plotted on a to 1 scale, 
with 1 representing the present values (which tantamounts to cu- 
mulating the equivalent Xt and yuj vectors for each region). The 
nucleus reached 80% of its mass at ^ = 8.5 Gyr, while the 7? > 1 
HLR region did so later, at ^ = L9 Gyr. As shown by Perez et 
al. (2013), this inside-out ordering of the mass assembly history 
applies to essentially all massive galaxies. 

Fig. [TT^ shows STUif). The SFR per unit area decreases 
from the nucleus outwards through most of the time, but the 
trend is reversed in the last ~ 500 Myr. Notice that, because 
of the smoothing, STUif) is now a continuous function of time, 
so that this panel represent the instantaneous SFR, as opposed to 



the mass-over-time definition in eq.[TT] Fig. [TT]f shows Scalos's 
b parameter for each region, which does use the running mean 
S'F'Rit) of eq. [llj This plot helps interpreting Fig. [9J which 
opens up the radial regions into full 2D maps, but compress the 
time axis by stipulating fixed tsF time scales. 



5.9. The whole versus the sum of the parts 

One of the applications of CALIFA is to assess the eff'ects of the 
lack of spatially resolved information on the derivation of phys- 
ical properties and SFHs from integrated- spectra surveys. The 
horizontal dashed magenta lines in Fig. 10 illustrate this point. 
They represent the results obtained from the analysis of the spec- 
trum of datacube as a whole, ie., adding upp all spaxels to em- 
ulate an integrated spectrum. Qualitatively, one expects that this 
should produce properties typical of the galaxy zones as a whole. 
This expectation is borne out in Fig.[TOJ where one sees that the 
extinction, mean ages and metallicities marked by the dashed 
magenta line do represent an overall average of the spatially re- 
solved values. In this particular example, they all match quite 
well the values at 1 HLR. The stellar masses obtained from the 
whole and the sum of the parts are also in excellent agreement, 
diff'erening by just 1%. Similarly, galaxy- wide luminosity and 
mass weighted mean ages and metallicities computed in these 
two ways agree to within 0.05 dex. 

But can one derive the SFH of a galaxy out of an integrated 
spectrum? Fig. 11 says that, at least in CALIFA 277, the answer 
is yes. The dashed magenta and the black solid lines are very 
similar in all panels of this figure. As before, the former repre- 
sents the results obtained from the analysis of the spatially com- 
pressed datacube, while the black curves are computed adding 
the STARLIGHT rcsults for each spaxel. We are thus comparing the 
whole against the sum of the parts, and they match. 

At first sight, this may seem a trivial result, but it is not. 
Formally, one only expects this to happen if Ay is the same at all 
positions, a condition which is approximately meet in CALIFA 
277. In objects where dust has strong spatial variations, the spec- 
tral fitting of the global spectrum with a single Ay will inevitably 
operate compensations, like increasing the age to compensate 
for an understimated Ay and vice- versa. Also, the combined ef- 
fects of population gradients and kinematics (e.g., disc popula- 
tion contributing more to the wings of absorption lines than the 
ones from the bulge) are hard to predict. Since we are presenting 
results for a single example galaxy, it remains to be seen how 
general this "coincidence" is. 



^ For clarity, the cumulative mass is computed with the total mass 
which has participated in star formation, instead of the one corrected 
for returned mass. In other words, we cumulate M' instead of M. 



5.10. Space x time diagrams: SFHs in 2D 

As is evident at this point of the paper, one of the main chal- 
lenges involved in the analysis of the multi dimensional data 
built from the combination of the spatial dimensions with the 
t and Z dimensions opened by population synthesis is how to vi- 
sualize the results. All examples shown so far project (or average 
over) two or more of these axes. 

Fig. [12] shows an attempt to visualize galaxy evolution as 
a function of both space and time. The trick is to compress xy 
into R, and collapse the Z axis, producing a radially averaged 
SFH(7?, t) map. The left panel shows the luminosity density X 
at each radial position R and for each age t. As for the other 
panels in this figure, the original Jltzyx array from which this 
map was derived was smoothed in log t, marginalized over Z and 
collapsed into one spatial dimension. Unlike in Fig. 10 we now 
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Fig. 11. Time evolution of light and mass representations of the SFH, as derived from the starlight fits. The discrete Xtzxy and yuj^ 
"population vectors" involved in these curves were smoothed in log t, marginalized in Z and integrated over the spatial regions 
indicated. Dotted grey lines indicate the nucleus (central pixel). Green, red, and blue lines correspond to R < 0.5 HLR, R < I HLR 
and R > I HLR regions, respectively. The black line represents the whole galaxy evolution as reconstructed from the sum of its 
parts, while (as in Fig.fTOl) the magenta dashed line is used to represent the results obtained for the spatially integrated spectrum. 



use the area averaging method discussed in Section 5.7 but the 
results do not depend strongly on this choice. 

The £R^t diagram shows that the light from old stars is con- 
centrated in the bulge, while the youngest stars shine ~ equally 
at all R. Also, t ^ I Gyr stars are more smoothly distributed in 
ages in the disc than in the bulge, where they are concentrated 
in two populations, one at ~ 2.5 Gyr and the other older than 
10 Gyr. Fig. [12]) shows the formed stellar mass surface density 
as a function of R and t. As usual, due to the highly non-linear 
relation between light and mass in stars, the young populations 
which show up so well in light practically disappear when seen 
in mass. Young stars reappear in the right panel, which shows 
the evolution in time and space of the SFR surface density. As in 
Fig.pJ^, these are obtained from the instantaneous SFR at each 
spaxel, computed as in Asari et al. (2007): 



SFR(0 



m; _ logg m; 

A^ ~ Alog^ t 



(15) 



The age sampling is logarithmic (A log t = constant), so that the 
SFR is proportional to the stellar mass formed at t divided by t. 



The ST^R^t and AVr^i diagrams therefore carry the same infor- 
mation; it is the t~^ factor which makes them look so diff'erent. 

When interpreting these or any other plot involving age, one 
should always keep in mind the logarithmic age resolution. For 
instance, taken at face value, the spatially integrated SFR of 
CALIFA 277 in the last few Myr is almost equal to its histor- 
ical peak, ~ 2.5 Gyr ago (solid line in Fig. 12:). However, these 
two peaks span vastly different time intervals, roughly by a fac- 
tor of Gyr/Myr = a thousand. The peak SFR around 2.5 Gyr ago, 
and indeed even well before that (at > 10 Gyr, when most of 
the mass was formed), were surely higher than the one we are 
now seeing at a few Myr, but cannot be resolved in time. No fos- 
sil record method will ever be able to distinguish bursts much 
shorter than their current age. These limitations are well known 
in the field, but it is fit to recall them to avoid missinterpretation 
of the results. 



Notwithstanding such age resolution issues. Fig. [12] repre- 
sents a novel way of looking at galaxies both in time and space, 
and it shows clear evidence in favor of an inside out growth 
scenario, in which the outer regions assemble their mass at a 
slower pace and a over more extended period. This behavior is 
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Fig. 12. R-t diagrams showing the radially averaged distribution of light, mass and SFR as a function distance from the nucleus and 
age. Left: Luminosity at /I = 5635 A per unit area. (Jlu^t)- Middle: Stellar mass formed per unit area (Al^ ^). Right: Time dependent 
SFR per unit area (ST'R). The solid (white) lines represent the sum over all spaxels for a given age. 
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Fig. 13. "Pseudo snapshots" of the cumulative stellar mass sur- 
face density of CALIFA 277 at different ages: From logt [yr] 
= 10.15 to 8 (as labeled in each panel). 



not unique to this galaxy. Perez et al (2012) find just the same in 
a study of the first 105 galaxies observed by CALIFA. 



5.11. Space vs. time "snapshots" 

There is so much one can plot in a single image. To go beyond 
the R-t diagrams of Fig. [12] one needs more dimensions than 
a sheet of paper can accommodate. An alternative is shown in 



Fig. 13 which shows ^slices through the cumulative M'^yf cube. 
UnliEe the previous diagrams, which compressed the SFH in one 
way or another, these plots reveal the full 2D nature of the mass 
assembly history of CALIFA 277. 



Age sequences like those in Fig. [13] can be constructed for 
several, but not all properties inferred though starlight or any 
other fossil method. One can, for instance, replace mass, by in- 
trinsic luminosity, SFR or stellar metallicity, in cumulative or 
diff'erential form, but it is obviously impossible to reconstruct 
maps of Ay, cr^, v* as a function of t. It is also worth pointing 
out that the panels in Fig. 13 are not truly snapshots of a movie, 
ie., they are not pictures of the galaxy as it appeared at diff'er- 
ent look-back times. Instead, these are maps of where stars of a 
given age t are located nowadays. 

Simulators should take notice of these simple facts. IFS data 
plus fossil methods provide a rich, yet inevitably limited form 
of time-travel. Illustrative and beautiful as they are, movies of 
stars and gas particles moving through time and space will never 
be directly compared to anything observational. In other words, 
simulations should be convolved through this "reality filter". The 
observationally relevant predictions are the distribution of stellar 
ages and metallicities as a function of xy. 

5. 12. Emission line maps 

Independent of the stellar population applications described 
throughout this paper, the starlight fits provide a residual spec- 
trum Rx = Ox- Mx which is, inasmuch as possible, free of stel- 
lar light. This is of great aid in the measurement of emission 
lines, and hence in the derivation of a series of diagnostics of the 
warm ISM properties, such as its kinematics, nebular extinction 
and chemical abundance. The complementarity of nebular and 
stellar analysis has been amply explored in studies of integrated 
spectra (eg, Cid Fernandes et al. 2007; Asari et al. 2007; Perez- 
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Fig. 14. Residual spectral cube in a 3D rendering. 



Montero et al. 2010; Garcia-Benito & Perez-Montero 2012), and 
the addition of two spatial dimensions to this analysis holds great 
promise. Kehrig et al. (2012), for instance, combine starlight 
analysis and nebular emission to study the impact of post-AGB 
populations on the ionization of the ISM in early type galaxies. 
Fig. 14 shows a 7?^,xy residuals cube, with some of the main 
emissionnnes marked. The plot shows that emission lines are 
brighter towards the outside, in particular over the spiral arms, 
as could be guessed from the distribution of actively star-forming 
regions traced by our SFH analysis (e.g.. Fig. [9]). Despite its ac- 
tive nucleus, the nuclear line fluxes are not particularly bright, 
implying a low luminosity AGN. 



6. Summary 

This paper described a complete "analysis pipeline" which pro- 
cesses reduced IFS data through the starlight spectral synthe- 
sis code. Data for the spiral galaxy NGC 2916 (CALIFA 277) 
were used to illustrate the journey from (x,y,A) datacubes to 
multi-dimensional maps of the synthesis products, focusing on 
the technical and methodological aspects. The three main steps 
in the analysis can be summarized as follows: 

1. Pre-processing: This step comprises (a) construction of a 
spatial mask; (b) adjustments in the flag and error spectra; 
(c) rest- framing and resampling; (d) spatial binning. The last 
three stages were coded in a generic python package (qbick) 
which can run in a fully automated fashion. Spatial binning 
was accomplished with an Voronoi tesselation scheme, tuned 
to produce zones with S /N >20 spectra. In practice, because 
of the good quality of the data, no binning was necessary 
within R ^ I HLR. The issue of correlated errors in spatial 
binning was discussed and an empirical recipe to correct for 
it was presented. The output of these pre-processing steps are 
datacubes with fluxes, errors and flags ready for a detailed A- 
hy-A spectral analysis, as well as spatial masks, a zone map 
and other ancillary products. 

2. STARLIGHT fits and organization of the synthesis results: 
Spectral fits for all spatial zones were performed using a base 
of SSP spectra from a combination of MILES and Granada 
models, spanning ages from 1 Myr to 14 Gyr and 4 metallic- 
ities between 0.2 and 1.5 solar. The results were packed in a 
coherent multi-layered FITS (or HDF5) file with the Python 
Calif a Starlight Synthesis Organizer, pycAsso. 

3. Post-processing tools: pycAsso includes a reader module and 
a series of analysis tools to perform operations like zone-to- 



pixel image reconstruction, mapping the spectral and stel- 
lar population properties derived by starlight into multi- 
dimensional time, metallicity, and spatial coordinates, aver- 
aging in spatial and temporal dimensions, manipulation of 
the stellar population arrays, etc. 

Some of the products of this whole analysis are standard in 
IFS work, like the kinematical field and emission line maps. The 
real novelty, of course, resides in the spatially resolved stellar 
population products recovered from the fossil record of galaxy 
evolution encoded in the datacubes. The 2D products range from 
maps of the stellar extinction and surface densities of stellar 
mass to more elaborate products like mean ages and metallici- 
ties, time averaged star formation rates and /^-parameters. These 
quantities are all normally used in the analysis of integrated 
spectra, but their application to IFS data raises some conceptual 
issues, related to the unavoidable fact that not all stars in a given 
spaxel were born there. On the other hand, one can use spatially 
resolved data to construct diagnostics unaplicable to integrated 
data, such as indices comparing the present and past SFR in dif- 
ferent regions, or comparing the "present here" with the "present 
elsewhere". These IFS -based variations of the traditional b pa- 
rameter were shown to be useful to highlight diff'erent aspects of 
the spatially resolved SFH. 

ID profiles in temporal or radial coordinates were used as 
a means to help interpreting the results. In our example galaxy, 
these compressed views of the (t, Z, x, y) manifold reveal clear 
age and metallicity negative gradients, as well as spatially de- 
pendent mass growth speeds compatible with an inside-out for- 
mation scenario. Finally, radius-age diagrams and "snapshot" se- 
quences were introduced as further means of visualization. 

Clearly, the application of spectral synthesis methods to IFS 
datacubes off'ers new ways of studying galaxies and their evo- 
lution. Future communications will use the tools laid out in this 
paper to explore the astrophysical implications of the results for 
galaxies in the CALIFA survey. 
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